#include "Point.h"
#include "Geometry.h"
#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <vector>

#define TEMPLATE template <size_t DIM, typename T>

TEMPLATE class Mesh;

TEMPLATE class IrregularMesh;

class Easymesh;
class RegularMesh;

TEMPLATE class Mesh
{
public:
    virtual size_t get_n_points() const = 0;
    virtual size_t get_n_edges() const = 0;
    virtual size_t get_n_grids() const = 0;
    virtual size_t get_n_dofs() const = 0;
    virtual const Point<DIM, T> &get_point(size_t _idx) const = 0;
    virtual const Geometry &get_edge(size_t _idx) const = 0;
    virtual const Geometry &get_grid(size_t _idx) const = 0;
    virtual std::vector<int> get_all_boundary() { };
    virtual Point<2, double> &get_dof(size_t _idx) = 0;
    virtual bool IsBoundary(size_t _idx){};
    //下面的都是给结构网格的接口，非结构网格不可用
    virtual const Point<2, double> &get_lbc() const {};
    virtual const Point<2, double> &get_ruc() const {};
    virtual size_t get_nx() const {};
    virtual size_t get_ny() const {};
    virtual size_t CortoIdx(int i, int j) {};
    virtual Point<DIM, T> &get_point(size_t _idx) =0;
    virtual Geometry &get_edge(size_t _idx) =0;
    virtual Geometry &get_grid(size_t _idx) =0;
};

class Easymesh : public Mesh<2, double>
{
private:
    std::valarray<Point<2, double> > nodes;
    std::valarray<Geometry> edges;
    std::valarray<Geometry> grids;
    std::valarray<Point<2, double> > dofs;
public:
    void readData(const std::string &filename);
    Easymesh() : Mesh() {};
    Easymesh(const std::string &filename);
    size_t get_n_points() const;
    size_t get_n_edges() const;
    size_t get_n_grids() const;
    size_t get_n_dofs() const;
    const Point<2, double> &get_point(size_t _idx) const;
    const Point<2, double> &get_dof(size_t _idx) const;
    const std::valarray<Point<2, double> > &get_dofs() const;
    const Geometry &get_edge(size_t _idx) const;
    const Geometry &get_grid(size_t _idx) const;
    Point<2, double> &get_point(size_t _idx);
    Point<2, double> &get_dof(size_t _idx);
    std::valarray<Point<2,double> > &get_dofs();
    Geometry &get_edge(size_t _idx);
    Geometry &get_grid(size_t _idx);
    void build_dofs();
    bool IsBoundary(size_t _idx);
};

const Point<2, double> &Easymesh::get_point(size_t _idx) const
{
    return nodes[_idx];
};
Point<2, double> &Easymesh::get_point(size_t _idx)
{
    return nodes[_idx];
};
const Point<2, double> &Easymesh::get_dof(size_t _idx) const
{
    return dofs[_idx];
};
Point<2,double> &Easymesh::get_dof(size_t _idx)
{
    return dofs[_idx];
};
const std::valarray<Point<2, double> > &Easymesh::get_dofs() const
{
    return dofs;
};
std::valarray<Point<2, double> > &Easymesh::get_dofs()
{
    return dofs;
};
const Geometry &Easymesh::get_edge(size_t _idx) const
{
    return edges[_idx];
};
Geometry &Easymesh::get_edge(size_t _idx)
{
    return edges[_idx];
};
const Geometry &Easymesh::get_grid(size_t _idx) const
{
    return grids[_idx];
};
Geometry &Easymesh::get_grid(size_t _idx)
{
    return grids[_idx];
};
size_t Easymesh::get_n_points() const
{
    return nodes.size();
};
size_t Easymesh::get_n_dofs() const
{
    return dofs.size();
};
size_t Easymesh::get_n_edges() const
{
    return edges.size();
};
size_t Easymesh::get_n_grids() const
{
    return grids.size();
};
bool Easymesh::IsBoundary(size_t _idx)
{
    return (bool) dofs[_idx].get_boundary_mark();
};
void Easymesh::readData(const std::string &filename)
{
    std::cout<< "Reading easymesh data file ..." <<std::endl;
    std::cout<< "\t reading node data ... " << std::flush;//强行输出完数据
    std::ifstream is((filename + ".n").c_str());
    size_t n_node, n_side, n_element;
    char text[64];
    is >> n_node >> n_element >> n_side;
    is.getline(text, 64);
    nodes.resize(n_node);
    for(size_t i=0; i < n_node; i++)
    {
        size_t ivtx;
        Geometry::bm_t bm;
        is >> ivtx >> nodes[i][0] >> nodes[i][1] >>bm;
        nodes[i].set_index(ivtx);
        nodes[i].set_boundary_mark(bm);
        nodes[i].set_vertex(0, ivtx);
        nodes[i].set_boundary(0, ivtx);
        nodes[i].set_dof(0, ivtx);
    }
    is.close();
    dofs = nodes;
    std::cout<< "OK!" <<std::endl;
    std::cout<< "\t reading side data ..." <<std::flush;
    is.open((filename + ".s").c_str());
    size_t n_edge;
    is >> n_edge;
    if(n_edge != n_side)
    {
        std::cerr <<"in side file: side number error. "<< std::endl;
        exit(-1);
    }
    edges.resize(n_side);
    for(size_t i=0; i<n_side ;i++)
    {
        size_t iedg, a, b;
        int ea, eb;
        Geometry::bm_t bm;
        edges[i].set_n_vertex(2);
        edges[i].set_n_boundary(2);
        edges[i].set_n_neighbor(2);
        edges[i].set_n_dofs(2);
        is >> iedg >> a >> b >> ea >> eb >> bm;
        edges[i].set_index(iedg);
        edges[i].set_boundary_mark(bm);
        edges[i].set_vertex(0, a);
        edges[i].set_vertex(1, b);
        edges[i].set_boundary(0, a);
        edges[i].set_boundary(1, b);
        edges[i].set_neighbor(0, ea);
        edges[i].set_neighbor(1, eb);
        edges[i].set_dof(0, a);
        edges[i].set_dof(1, b);
    }
    is.close();
    std::cout<< "OK!" <<std::endl;
    std::cout<<"\t reading element data ..." <<std::flush;
    is.open((filename + ".e").c_str());
    size_t n_ele;
    is >> n_ele;
    if(n_ele != n_element)
    {
        std::cerr<< "in element file : element number error."<< std::endl;
        exit(-1);
    } 
    size_t n_n;
    is >> n_n;
    if(n_n != n_node)
    {
        std::cerr <<" in element file: node number error." << std::endl;
        exit(-1);
    }
    size_t n_s;
    is >> n_s;
    if(n_s != n_side)
    {
        std::cerr << "in element file : side number error." << std::endl;
        exit(-1);
    }
    is.getline(text, 64);
    grids.resize(n_element);
    for(size_t i=0; i< n_element; i++)
    {
        grids[i].set_n_vertex(3);
        grids[i].set_n_boundary(3);
        grids[i].set_n_neighbor(3);
        size_t igrd, v0, v1, v2, b0, b1, b2;
        int n0, n1, n2;
        is >> igrd >> v0 >> v1 >> v2 >>n0 >> n1 >> n2 >>b0 >>b1 >> b2;
        grids[i].set_index(igrd);
        grids[i].set_boundary_mark(0);
        grids[i].set_vertex(0, v0);
        grids[i].set_vertex(1, v1);
        grids[i].set_vertex(2, v2);
        grids[i].set_neighbor(0, n0);
        grids[i].set_neighbor(1, n1);
        grids[i].set_neighbor(2, n2);
        grids[i].set_boundary(0, b0);
        grids[i].set_boundary(1, b1);
        grids[i].set_boundary(2, b2);
        const Point<2, double> &p0 = nodes[v0];
        const Point<2, double> &p1 = nodes[v1];
        const Point<2, double> &p2 = nodes[v2];
        if((p1[0] - p0[0])*(p2[1]-p0[1])-(p1[1] -p0[1])*(p2[0]-p0[0]) < 0)//三角形的面积
        {
            std::cerr << "vertices of element "<< igrd << "is not correctly ordered. "<< std::endl;
            exit(-1);
        }
        grids[i].set_n_dofs(3);
        grids[i].set_dof(0, v0);
        grids[i].set_dof(1, v1);
        grids[i].set_dof(2, v2);
    }
    is.close();
    std::cout << " OK! "<< std::endl;
};

Easymesh::Easymesh(const std::string &filename) : Mesh()
{
    readData(filename);
};

void Easymesh::build_dofs()
{
    int dof_index = 0;
    int n_ele = get_n_grids();
    int n_total_dofs = get_n_points() + get_n_edges();
    dofs.resize(n_total_dofs);
    for(int i =0; i< n_ele; i++)
    {
        int n_local_dofs = 6;
        grids[i].set_n_dofs(n_local_dofs);
        int n_vtx = grids[i].get_n_vertex();
        for(int k =0; k < n_vtx; k++)
        {
            Point<2, double> &vtx = get_point(grids[i].get_vertex(k));
            vtx.set_dof(0, -1);
        }
        int n_bnd = grids[i].get_n_boundary();
        for(int k =0; k < n_bnd; k++)
        {
            Geometry &bnd = get_edge(grids[i].get_boundary(k));
            if(bnd.get_n_dofs() == 2)
            {
                int a =bnd.get_dofs(0);
                int b =bnd.get_dofs(1);
                bnd.set_n_dofs(3);
                bnd.set_dof(0, a);
                bnd.set_dof(1, b);
                bnd.set_dof(2, -1);
            }
        }
    }
    for(int i=0; i < n_ele; i++)
    {
        int n_local_dofs = grids[i].get_n_dofs();
        int n_vtx = grids[i].get_n_vertex();
        for(int k =0; k < n_vtx; k++)
        {
            Point<2, double> &vtx = get_point(grids[i].get_vertex(k));
            int global_dof = vtx.get_dofs(0);
            if(global_dof == -1)
            {
                vtx.set_dof(0, dof_index);
                grids[i].set_dof(k, dof_index);
                vtx.set_index(dof_index);
                dofs[dof_index] = vtx;
                dof_index ++;
            }
            else
                grids[i].set_dof(k, global_dof);
        }
        int n_bnd = grids[i].get_n_boundary();
        for(int k =0; k < n_bnd; k++)
        {
            Geometry &edge = get_edge(grids[i].get_boundary(k));
            int global_dof = edge.get_dofs(2);
            if(global_dof == -1)
            {
                edge.set_dof(2, dof_index);
                grids[i].set_dof(k + n_vtx, dof_index);
                Point<2, double> &vtx0 = get_point(edge.get_dofs(0));
                Point<2, double> &vtx1 = get_point(edge.get_dofs(1));
                Point<2, double> vtx({(vtx0[0]+vtx1[0])/2, (vtx0[1]+vtx1[1])/2});
                vtx.set_dof(0, dof_index);
                vtx.set_boundary_mark(edge.get_boundary_mark());
                vtx.set_index(dof_index);
                dofs[dof_index] = vtx;
                dof_index ++;
            }
            else
                grids[i].set_dof(k + n_vtx, global_dof);
        }
    }
};

class RegularMesh : public Mesh<2, double>
{
protected:
    Point<2, double> lbc; // x0y0
    Point<2, double> ruc; // x1y1
    size_t nx;
    size_t ny;

public:
    Point<2, double> node;
    Geometry edge;
    Geometry grid;

    void Initial();
    RegularMesh();
    RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny);
    void readData(const std::string &filename){ };//do nothing
    const Point<2, double> &get_lbc() const { return lbc; };
    const Point<2, double> &get_ruc() const { return ruc; };
    void set_lbc(const Point<2, double> &_lbc){ lbc = _lbc; };
    void set_ruc(const Point<2, double> &_ruc){ ruc = _ruc; };
    size_t get_nx() const { return nx; };
    size_t get_ny() const { return ny; };
    double get_hx() const { return (ruc[0] - lbc[0]) / nx ;};
    double get_hy() const { return (ruc[1] - lbc[1]) / ny ;};
    void set_nx(size_t _nx) { nx = _nx; };
    void set_ny(size_t _ny) { ny = _ny; };
    void set_nxny(size_t _nx, size_t _ny) { nx = _nx; ny = _ny; };
    size_t get_n_points() const;
    size_t get_n_edges() const;
    size_t get_n_grids() const;
    size_t get_n_dofs() const;
    Point<2, double> &get_point(size_t _idx);
    virtual Point<2, double> &get_dof(size_t _idx);
    Geometry &get_edge(size_t _idx);
    virtual Geometry &get_grid(size_t _idx);
    void _get_grid(size_t _idx);
    std::vector<int> get_all_boundary();
    bool IsBoundary(size_t _idx)
    {
        std::vector<int> _BndIndex = this->get_all_boundary();
        return !(std::find(_BndIndex.begin(), _BndIndex.end(), _idx) == _BndIndex.end());
    };
    virtual size_t CortoIdx(int i, int j){ return (nx + 1) * i + j ;};
            //以下缺少相关实现
    virtual const Point<2, double> &get_point(size_t _idx) const { };
    virtual const Geometry &get_edge(size_t _idx) const {};
    virtual const Geometry &get_grid(size_t _idx) const {};
};

void RegularMesh::Initial()
{
    edge.set_n_dofs(2);
    edge.set_n_vertex(2);
    edge.set_n_neighbor(2);
    grid.set_n_neighbor(3);
    grid.set_n_vertex(3);
    grid.set_n_dofs(3);
};
RegularMesh::RegularMesh() : Mesh()
{
    lbc[0] = lbc[1] = 0;
    ruc[0] = ruc[1] = 1;
    nx = ny = 2;
    Initial();
};

RegularMesh::RegularMesh(const Point<2, double> &_lbc, const Point<2, double> &_ruc, int _nx, int _ny) : Mesh(),
            lbc(_lbc), ruc(_ruc), nx(_nx), ny(_ny) { Initial(); }

Point<2, double>  &RegularMesh::get_dof(size_t _idx){ return get_point(_idx); };

Geometry &RegularMesh::get_edge(size_t _idx)
{
    size_t ix = (_idx) % (3* nx +1);
    size_t iy = (_idx) / (3 *nx +1);
    edge.set_index(_idx);
    if(ix < nx)
    {
        edge.set_neighbor(0, ((iy == 0)? -1:(iy - 1)*(2 *nx) +2 *ix));
        edge.set_neighbor(1, ((ix == ny)? -1: iy*2*nx + 2 * ix +1 ));
        edge.set_dof(0, iy * (nx +1) + ix);
        edge.set_dof(1, iy *(nx +1)+ ix +1);
        edge.set_vertex(0, iy * (nx +1) + ix);
        edge.set_vertex(1, iy *(nx +1) + ix +1);
    }
    else
    {
        ix -= nx;
        if(ix %2 == 0)
        {
            edge.set_neighbor(0, ((ix == 0) ? -1 : iy * (2 * nx ) +2 * ix -1));
            edge.set_neighbor(1, ((ix > nx) ? -1 : iy * (2 * nx) + 2 * ix));
            edge.set_dof(0, iy * (nx +1) + ix);
            edge.set_dof(1, (iy +1) * (nx +1) + ix);
            edge.set_vertex(0, iy * (nx +1) + ix);
            edge.set_vertex(1, (iy +1) * (nx +1) + ix);
        }
        else
        {
            ix /= 2;
            edge.set_neighbor(0, iy *(2 * nx) + 2 * ix);
            edge.set_neighbor(1, iy *(2 *nx) + 2 *ix +1 );
            edge.set_dof(0, iy *(nx +1) + ix +1);
            edge.set_dof(1, (iy +1) * ( nx +1) + ix);
            edge.set_vertex(0, iy *( nx +1) + ix +1);
            edge.set_vertex(1, (iy +1) * (nx +1) + ix);
        }
    }
    edge.set_boundary_mark((ix == 0 || iy == 0 || ix == nx || iy== ny));
    return edge;
};

Geometry &RegularMesh::get_grid(size_t _idx)
{
    _get_grid(_idx);
    return grid;
};

void RegularMesh::_get_grid(size_t _idx)
{
    size_t dof[3];
    size_t en[3];
    grid.set_index(_idx);
    if(_idx % 2 == 0)
    {
        size_t ix = (_idx / 2) % nx;
        size_t iy = (_idx /2) / nx;
        dof[0] = iy * (nx +1) +ix;
        dof[1] = iy * (nx +1) + ix +1;
        dof[2] = (iy +1) * (nx +1) +ix;
        en[0] = _idx +1;
        en[1] = (ix == 0) ? -1:_idx -1;
        en[2] = (iy == 0) ? -1:_idx -2 * nx +1;
    }
    else
    {
        size_t ix = ((_idx -1) / 2) % nx;
        size_t iy = ((_idx -1) / 2) / nx;
        dof[0] = iy * (nx +1) +ix +1;
        dof[1] = (iy + 1) *(nx +1) +ix +1;
        dof[2] = (iy +1)*(nx +1) +ix;
        en[0] = (iy == ny)? -1:_idx +2*nx-1;
        en[1] = _idx -1;
        en[2] = (ix == nx)? -1:_idx +1;
    }
    for(int i = 0; i< 3; i++) 
    {
        grid.set_vertex(i, dof[i]);
        grid.set_dof(i, dof[i]);
        grid.set_neighbor(i, en[i]);
    }
};
size_t RegularMesh:: get_n_points() const { return (nx +1) *(ny +1);}
size_t RegularMesh::get_n_dofs() const { return (nx +1)*(ny+1);}
size_t RegularMesh::get_n_edges() const { return (nx +1) * ny + (ny +1)* nx +nx*ny ;}
size_t RegularMesh::get_n_grids() const { return nx * ny *2 ;}
std::vector<int> RegularMesh::get_all_boundary()
{
    std::vector<int> _BndIndex;
    for(int i = 0; i <= nx; i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs() - i - 1);
    }
    for(int j = 0; j <= ny -1; j++)
    {
        _BndIndex.push_back(j *(nx +1));
        _BndIndex.push_back((j+1) * (nx +1) -1);
    }
    return _BndIndex;
};

class RegularMeshP2 : public RegularMesh
{
public:
    RegularMeshP2() : RegularMesh() { grid.set_n_dofs(6);};
    RegularMeshP2(const Point<2, double> & _lbc, const Point<2, double> &_ruc, int _nx, int _ny):
        RegularMesh(_lbc, _ruc, _nx, _ny){ grid.set_n_dofs(6); };
    Point<2, double> &get_dof(size_t _idx);
    Geometry &get_grid(size_t _idx);
    size_t get_n_dofs() const;
    std::vector<int> get_all_boundary();
    size_t CortoIdx(int i, int j) { return (2 * nx +1)*i + j ;};
};

Point<2,double> &RegularMeshP2::get_dof(size_t _idx)
{
    size_t ix = (_idx) % ( 2 *nx +1);
    size_t iy = (_idx) /(2 * nx +1);
    double hx = get_hx(), hy = get_hy();
    node[0] = lbc[0] + ix * hx /2.0;
    node[1] = lbc[1] + iy * hy /2.0;
    node.set_index(_idx);
    node.set_n_dofs(1);
    node.set_n_vertex(1);
    node.set_vertex(0, _idx);
    node.set_boundary_mark((ix == 0 || iy == 0 || ix == 2 *nx || iy == 2 *ny));
    return node;
};

Geometry &RegularMeshP2::get_grid(size_t _idx)
{
    size_t dof[6];
    _get_grid(_idx);
    if(_idx % 2 == 0)
    {
        size_t ix = (_idx / 2) % nx;
        size_t iy = (_idx / 2) / nx;
        dof[0] = 2 * iy * ( 2 * nx +1) + 2 * ix;
        dof[1] = 2 * iy *(2 * nx +1) + 2 *(ix +1);
        dof[2] = 2 * (iy +1) * (2 * nx +1) + 2 *ix;
        dof[3] = 2 * (iy + 0.5) * (2 * nx +1) +2 *(ix + 0.5);
        dof[4] = 2 * (iy + 0.5) * (2 * nx +1) + 2 * ix;
        dof[5] = 2 * iy * (2 * nx +1) +2 * (ix + 0.5);
    }
    else 
    {
        size_t ix = ((_idx - 1) / 2) % nx;
        size_t iy = ((_idx - 1) / 2) / nx;
        dof[0] = 2 * iy *(2 *nx +1) +2 *(ix +1);
        dof[1] = 2 *(iy +1)* (2 *nx +1) + 2 *(ix +1);
        dof[2] = 2 *(iy + 1) * (2 * nx +1) + 2 * ix;
        dof[3] = 2 *(iy +1) * (2 *nx +1) + 2*(ix + 0.5);
        dof[4] = 2 *(iy +0.5) *(2 * nx +1) + 2 *(ix +0.5);
        dof[5] = 2 * (iy + 0.5) * (2 *nx +1) + 2 *(ix +1);
    }
    for(int i = 0; i < 6; i++)
        grid.set_dof(i, dof[i]);
    return grid;
}

size_t RegularMeshP2:: get_n_dofs() const
{
    return (2 *nx +1) *(2 *ny +1);
}
std::vector<int> RegularMeshP2::get_all_boundary()
{
    std::vector<int> _BndIndex;
    for(int i =0; i<= 2 *nx; i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs() - i -1);
    }
    for(int j = 1; j <=2*ny -1; j++)
    {
        _BndIndex.push_back(j * (2 *nx +1));
        _BndIndex.push_back((j+1)*(2 *nx +1) -1);
    }
    return _BndIndex;
};

class RegularMeshP3 : public RegularMesh
{
public:
    RegularMeshP3() : RegularMesh() { grid.set_n_dofs(10);};
    RegularMeshP3(const Point<2, double> & _lbc, const Point<2, double> &_ruc, int _nx, int _ny):
        RegularMesh(_lbc, _ruc, _nx, _ny){ grid.set_n_dofs(10); };
    Point<2, double> &get_dof(size_t _idx);
    Geometry &get_grid(size_t _idx);
    size_t get_n_dofs() const;
    std::vector<int> get_all_boundary();
    size_t CortoIdx(int i, int j) { return (3 * nx +1)*i + j ;};
};

Point<2,double> &RegularMeshP3::get_dof(size_t _idx)
{
    size_t ix = (_idx) % ( 3 *nx +1);
    size_t iy = (_idx) / ( 3* nx +1);
    double hx = get_hx(), hy = get_hy();
    node[0] = lbc[0] + ix * hx /3.0;
    node[1] = lbc[1] + iy * hy /3.0;
    node.set_index(_idx);
    node.set_n_dofs(1);
    node.set_n_vertex(1);
    node.set_vertex(0, _idx);
    node.set_boundary_mark((ix == 0 || iy == 0 || ix == 3 *nx || iy == 3 *ny));
    return node;
};

Geometry &RegularMeshP3::get_grid(size_t _idx)
{
    size_t dof[10];
    _get_grid(_idx);
    if(_idx % 2 == 0)
    {
        size_t ix = (_idx / 2) % nx;
        size_t iy = (_idx / 2) / nx;
        dof[0] = 3 * iy * ( 3 * nx +1) + 3 * ix;
        dof[1] = 3 * iy *(3 * nx +1) + 3 *(ix +1);
        dof[2] = 3 * (iy +1) * (3 * nx +1) + 3 *ix;
        dof[3] = 3 * iy  * (3 * nx +1) +3 *ix + 1;
        dof[4] = 3 * iy  * (3 * nx +1) +3 *ix + 2;
        dof[5] = 3 * iy  * (3 * nx +1) +3 *ix + 2 + 3 *nx +1;
        dof[6] = (3 *iy + 2) * ( 3 * nx + 1) + 3 * ix + 1;
        dof[7] = (3 *iy + 2) * ( 3 * nx + 1) + 3 * ix ;
        dof[8] = (3 *iy + 1) * ( 3 * nx + 1) + 3 * ix ;
        dof[9] = (3 * iy + 1) * ( 3 * nx + 1) + 3 * ix +1;
    }
    else 
    {
        size_t ix = ((_idx - 1) / 2) % nx;
        size_t iy = ((_idx - 1) / 2) / nx;
        dof[0] = 3 * iy *(3 *nx +1) +3 *(ix +1);
        dof[1] = 3 *(iy +1)* (3 *nx +1) + 3 *(ix +1);
        dof[2] = 3 *(iy + 1) * (3 * nx +1) + 3 * ix;
        dof[3] = (3 * iy + 1) *(3 *nx +1) +3 *(ix +1);
        dof[4] = (3 * iy + 2) *(3 *nx +1) +3 *(ix +1);
        dof[5] = 3 *(iy + 1) *(3 * nx + 1) + 3 *ix +2;
        dof[6] = 3 *(iy + 1) *(3 * nx + 1) + 3 *ix +1;
        dof[7] = (3 *iy + 2) *(3 * nx + 1) + 3 *ix +1;
        dof[8] = (3 *iy + 1) *(3 * nx + 1) + 3 *ix +2;
        dof[9] = (3 *iy + 2) * (3 *nx +1) + 3 * ix + 2;
    }
    for(int i = 0; i < 10; i++)
        grid.set_dof(i, dof[i]);
    return grid;
}

size_t RegularMeshP3:: get_n_dofs() const
{
    return (3 *nx +1) *(3 *ny +1);
}
std::vector<int> RegularMeshP3::get_all_boundary()
{
    std::vector<int> _BndIndex;
    for(int i =0; i<= 3 *nx; i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs() - i -1);
    }
    for(int j = 1; j <=3*ny -1; j++)
    {
        _BndIndex.push_back(j * (3 *nx +1));
        _BndIndex.push_back((j+1)*(3 *nx +1) -1);
    }
    return _BndIndex;
};

class RegularMeshP4 : public RegularMesh
{
public:
    RegularMeshP4() : RegularMesh() { grid.set_n_dofs(15);};
    RegularMeshP4(const Point<2, double> & _lbc, const Point<2, double> &_ruc, int _nx, int _ny):
        RegularMesh(_lbc, _ruc, _nx, _ny){ grid.set_n_dofs(15); };
    Point<2, double> &get_dof(size_t _idx);
    Geometry &get_grid(size_t _idx);
    size_t get_n_dofs() const;
    std::vector<int> get_all_boundary();
    size_t CortoIdx(int i, int j) { return (4 * nx +1)*i + j ;};
};

Point<2,double> &RegularMeshP4::get_dof(size_t _idx)
{
    size_t ix = (_idx) % ( 4 *nx +1);
    size_t iy = (_idx) / ( 4* nx +1);
    double hx = get_hx(), hy = get_hy();
    node[0] = lbc[0] + ix * hx /4.0;
    node[1] = lbc[1] + iy * hy /4.0;
    node.set_index(_idx);
    node.set_n_dofs(1);
    node.set_n_vertex(1);
    node.set_vertex(0, _idx);
    node.set_boundary_mark((ix == 0 || iy == 0 || ix == 4 *nx || iy == 4 *ny));
    return node;
};

Geometry &RegularMeshP3::get_grid(size_t _idx)
{
    size_t dof[15];
    _get_grid(_idx);
    if(_idx % 2 == 0)
    {
        size_t ix = (_idx / 2) % nx;
        size_t iy = (_idx / 2) / nx;
        dof[0] = 4 * iy * ( 4 * nx +1) + 4 * ix;
        dof[1] = 4 * iy *(4 * nx +1) + 4 *(ix +1);
        dof[2] = 4 * (iy +1) * (4 * nx +1) + 4 *ix;
        dof[3] = (4 *iy +1)*(4 *nx +1)+ 4 *ix + 3;
        dof[4] = (4*iy +2) *(4 *nx +1) + 4 *ix +2;
        dof[5] = (4 *iy +3) *(4 *nx +1) + 4 *ix +1;
        dof[6] = (4 *iy +3) *(4 *nx +1) +4 *ix;
        dof[7] = (4 *iy +2) * (4 *nx +1) + 4 *ix;
        dof[8] = (4 *iy +1) * (4 *nx +1) + 4 *ix;
        dof[9] = 4 *iy *(4 *nx +1) +4 *ix +1;
        dof[10] = 4 *iy *(4 *nx +1) +4 *ix +2;
        dof[11] = 4 *iy *(4 *nx +1) +4 *ix +3;
        dof[12] = (4 *iy +2) *(4 *nx +1) + 4 *ix +1;
        dof[13] = (4 *iy +1) *(4 *nx +1) +4 *ix +1;
        dof[14] = (4 *iy +1) *(4 *nx +1) +4 *ix +2;
    }
    else 
    {
        size_t ix = ((_idx - 1) / 2) % nx;
        size_t iy = ((_idx - 1) / 2) / nx;
        dof[0] = 4 * iy *(4 *nx +1) +4 *(ix +1);
        dof[1] = 4 *(iy +1)* (4 *nx +1) + 4 *(ix +1);
        dof[2] = 4 *(iy + 1) * (4 * nx +1) + 4 * ix;
        dof[3] = 4 *(iy +1)* (4 *nx +1) + 4 *ix +3 ;
        dof[4] = 4 *(iy +1)* (4 *nx +1) + 4 *ix +2;
        dof[5] = 4 *(iy +1)* (4 *nx +1) + 4 *ix +1 ;
        dof[6] = (4 *iy +3)* (4 *nx +1) + 4 *ix +1 ;
        dof[7] = (4 *iy +2)* (4 *nx +1) + 4 *ix +2 ;
        dof[8] = (4 *iy + 1)*(4 *nx +1) +4* ix +3;
        dof[9] = (4 *iy + 1)*(4 *nx +1) +4* ix +4;
        dof[10] = (4 *iy + 2)*(4 *nx +1) +4* ix +4;
        dof[11] = (4 *iy + 3)*(4 *nx +1) +4* ix +4;
        dof[12] = (4 *iy +3) *(4 *nx +1) +4 *ix +2;
        dof[13] = (4 *iy +2) *(4 *nx +1) +4 *ix +3;
        dof[14] = (4 *iy + 3)*(4 *nx +1) +4* ix +3;
    }
    for(int i = 0; i < 15; i++)
        grid.set_dof(i, dof[i]);
    return grid;
}

size_t RegularMeshP3:: get_n_dofs() const
{
    return (4 *nx +1) *(4 *ny +1);
}
std::vector<int> RegularMeshP3::get_all_boundary()
{
    std::vector<int> _BndIndex;
    for(int i =0; i<= 4 *nx; i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(get_n_dofs() - i -1);
    }
    for(int j = 1; j <=4*ny -1; j++)
    {
        _BndIndex.push_back(j * (3 *nx +1));
        _BndIndex.push_back((j+1)*(3 *nx +1) -1);
    }
    return _BndIndex;
};

template<typename T>
class RectangleMesh : public Mesh<2, T>
{
private:
    std::valarray<Point<2, T> > nodes;
    std::valarray<Geometry> edges;
    std::valarray<Geometry> grids;
    std::valarray<Point<2, T> > dofs;

public:
    RectangleMesh() : Mesh<2, T>() {};
    RectangleMesh(const Point<2, T> &_lbc, const Point<2, T> &_ruc, int _nx, int _ny);
    virtual size_t get_n_points() const;
    virtual size_t get_n_edges() const;
    virtual size_t get_n_grids() const;
    virtual size_t get_n_dofs() const;
    virtual const Point<2, T> &get_point(size_t _idx) const;
    virtual const Point<2, T> &get_dof(size_t _idx) const;
    virtual const std::valarray<Point<2, T> > &get_dofs() const;
    virtual const Geometry &get_edge(size_t _idx) const;
    virtual const Geometry &get_grid(size_t _idx) const;
};

template <typename T>
const Point<2, T> &RectangleMesh<T>::get_point(size_t _idx) const
{
    return nodes[_idx];
};
template <typename T>
const Point<2, T> &RectangleMesh<T>::get_dof(size_t _idx) const
{
    return dofs[_idx];
};
template <typename T>
const std::valarray<Point<2, T> > &RectangleMesh<T>::get_dofs() const { return dofs; };

template <typename T>
const Geometry &RectangleMesh<T>:: get_edge(size_t _idx) const { return edges[_idx]; };

template <typename T>
const Geometry &RectangleMesh<T>::get_grid(size_t _idx) const { return grids[_idx]; };

template <typename T>
size_t RectangleMesh<T>::get_n_points() const { return nodes.size(); };

template <typename T>
size_t RectangleMesh<T>::get_n_dofs() const { return dofs.size(); };

template <typename T>
size_t RectangleMesh<T>::get_n_edges() const { return edges.size(); };

template <typename T>
size_t RectangleMesh<T>::get_n_grids() const { return grids.size(); };

template <typename T>
RectangleMesh<T>::RectangleMesh(const Point<2, T> &_lbc, const Point<2, T> &_ruc, int _nx, int _ny)
{
    T hx =(_ruc[0] - _lbc[0]) / _nx;
    T hy =(_ruc[1] - _ruc[0]) / _ny;
    nodes.resize((_nx +1) * (_ny +1));
    edges.resize(_nx * (_ny +1) + (_nx +1) * _ny + _nx * _ny); // 这里怎么画出来还是三角形网格？
    grids.resize(_nx * _ny *2);
    for(int ix =0; ix < _nx +1; ix++)
        for(int iy =0; iy < _ny +1; i++)
        {
            int node_idx = iy * (_nx +1) +ix;
            nodes[node_idx].set_index(node_idx);
            if(ix == 0 || ix == _nx || iy ==0 || iy == _ny)
                nodes[node_idx].set_boundary_mark(1);
            nodes[node_idx].set_vertex(0, node_idx);
            nodes[node_idx].set_boundary(0, node_idx);
            nodes[node_idx].set_dof(0, node_idx);
            nodes[node_idx][0] = ix * hx;
            nodes[node_idx][1] = iy * hy;
        }
    for(int ix = 0; ix < _nx; ix ++) //不太懂为什么矩形网格不是四个点和四条边组成的吗
        for(int iy = 0; iy < _ny; iy ++)
        {
            int even_ele_idx = (iy * _nx +ix) *2;
            grids[even_ele_idx].set_n_vertex(3);
            grids[even_ele_idx].set_n_boundary(3);
            grids[even_ele_idx].set_n_neighbor(3);

            grids[even_ele_idx].set_index(even_ele_idx);
            int v0 = iy * (_nx +1) + ix;
            int v1 = iy * (_nx +1) +ix +1;
            int v2 = (iy +1) * (_nx +1) +ix;
            grids[even_ele_idx].set_vertex(0, v0);
            grids[even_ele_idx].set_vertex(1, v1);
            grids[even_ele_idx].set_vertex(2, v2);
            grids[even_ele_idx].set_n_dofs(3);
            grids[even_ele_idx].set_dof(0, v0);
            grids[even_ele_idx].set_dof(1, v1);
            grids[even_ele_idx].set_dof(2, v2);

            int odd_ele_idx = (iy * _nx +ix) *2 +1;
            grids[odd_ele_idx].set_n_vertex(3);
            grids[odd_ele_idx].set_n_boundary(3);
            grids[odd_ele_idx].set_n_neighbor(3);

            grids[odd_ele_idx].set_index(odd_ele_idx);
            v0 = (iy + 1) *(_nx +1) +ix +1;
            v1 = (iy +1) * (_nx +1) + ix;
            v2 = iy * (_nx +1) + ix +1;
            grids[odd_ele_idx].set_vertex(0, v0);
            grids[odd_ele_idx].set_vertex(1, v1);
            grids[odd_ele_idx].set_vertex(2, v2);
            grids[odd_ele_idx].set_n_dofs(3);
            grids[odd_ele_idx].set_dof(0 , v0);
            grids[odd_ele_idx].set_dof(1 , v1);
            grids[odd_ele_idx].set_dof(2 , v2);
        }
    dofs = nodes;
};
